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t We study numerically the aging dynamics of the two-dimensional p-state clock model after a 

. quench from an infinite temperature to the ferromagnetic phase or to the Kosterlitz-Thouless phase. 

The system exhibits the general scaling behavior characteristic of non-disordered coarsening systems. 
CO ' For quenches to the ferromagnetic phase, the value of the dynamical exponents, suggests that the 

O^l , model belongs to the Ising-type universality class. Specifically, for the integrated response function 

x(t,s) — s~ ax f(t/s), we find a x consistent with the value a x — 0.28 found in the two-dimensional 
Ising model. 
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I. INTRODUCTION 



The phase-ordering kinetics of systems quenched from an high-temperature disordered state to an ordered phase 
or to a critical phase is characterized by the growth of a characteristic length L(t). In the late stage a power-law 
behavior 

g L(i)~<5 (1) 

■"O ■ sets in, where z is the dynamical exponent, and dynamical scaling [jj] is observed. Accordingly, configurations of the 
system at two subsequent times are statistically equivalent if lengths are measured in units of L(t), namely if the 
rescaled length x = r/L(t) is considered. This property is reflected by the analytical form of physical observables. 
In the quench to a critical point, for example, the equal-time order parameter correlation function obeys the scaling 
form 



where d is the spatial dimensionality and i] is the usual exponent of static critical phenomena. Scaling behaviors such 
as Eqs. (|1I2|I are generic for growth kinetics in non-frustrated, non-disordered systems and are observed regardless of 
different specific details. Concerning the values of the exponents, such as z or others entering different quantities, they 
\Q ' are expected to take the same value for systems belonging to the same non-equilibrium universality class. It is well 
known that systems undergoing a second order equilibrium phase transition can be classified into static equilibrium 
universality classes according to the value of their critical indices. These are found to depend only on a small set of 
parameters, such as space dimensionality or the number of components of the order parameter . In the same way, 
dynamic universality classes can also be introduced on the basis of the value of dynamic exponents. This subject 
has been thoroughly studied for the equilibrium critical dynamics 0, where the renormalization group provides the 
basic mechanism for scaling and universality, analogously to the static case. By contrast, the same subject is not well 
understood in far from equilibrium systems 0, E| • 
£j , In this paper we consider the phase-ordering kinetics of the clock model with a p-fold degenerate ground state 
I* ■ in two dimensions. For p < 4 this model has a single second order phase transition separating a disordered from 
.£h ! an ordered phase. For p > 5 a Kosterlitz-Thouless (KT) critical phase also exists and the phase transitions are of 
the Kosterlitz-Thouless type. We study numerically the model with p — 3 and p — 6 and a non-conserved order 
parameter, quenched to the ordered region and, for the case p = 6, also to the KT phase. In previous works 0, 
the growth law and the scaling of G(r, t) in quenches to the ordered region was analyzed; for quenches to a 
critical point, two time quantities were studied in Q- Here we extend these results presenting a global analysis of the 
scaling properties for quenches in the ordered and in the critical region, by considering one-time quantities, such as 
G(r, t), and two-time quantities, such as the autocorrelation function and the integrated response function. We find 
that the scaling forms expected for all these quantities in non-disordered systems are obeyed, although preasymptotic 
corrections are observed in the simulated range of times. We may then conclude that the clock model exhibits the 
generic scaling behavior characteristic of phase-ordering systems in the late stage of the dynamics. This calls for the 
question of the value of the dynamical exponents and the issue of their universality, namely whether the dynamics 
of the model is regulated by the same exponents of other coarsening systems and whether their value depends on p. 



2 



Since in the KT phase critical exponents depend continuously on temperature, this problem is pertinent to the quench 
to the ordered phase. Among the most usually considered two-dimensional statistical models of ferromagnetism, the 
Ising and the XY model are those to which the clock model can be naturally compared. These models are particular 
cases, with p = 2 and p — oo, of the clock model. The former has a scalar order parameter and a discrete symmetry. 
In the latter, the order parameter is a vector with N = 2 components and there is a continuous 0(2) symmetry. The 
symmetry group of the model, together with the spatial dimensionality, determines the types of topological defects. 
Consequently, topological defects are interfaces or two-dimensional vortices in the Ising and XY model, respectively. 
Since the nature of the topological defects controls several equilibrium properties and the late stage ordering |lj these 
models belong to different universality classes both in equilibrium and out of equilibrium. The clock model is, in some 
sense, intermediate between the Ising and XY model, because the order parameter is a vector with two components, 
as in the XY model, but there is a finite degeneracy of the ground state, namely a discrete symmetry. Topological 
defects are then both interfaces and vortices. 

In this paper we show that the exponents measured for the phase-ordering kinetics of the clock model with p = 3, 6 
are consistent with those of the two-dimensional Ising model quenched below the critical temperature. This suggests 
that systems with a finite degeneracy of the ground state may belong to the Ising non-equilibrium universality class, 
and that the presence of other topological defects besides interfaces does not affect the universal properties of these 
systems. 

This Article is organized as follows: In Sec. [H] we introduce the model and define the main observables. In Sec. IIIII 
we discuss the general scaling behavior of systems quenched into an ordered region, present the results of numerical 
simulations of the clock model with p — 3 and p — 6, and compare them with the behavior of the Ising model and of 
the XY model. In Sec. II VI we discuss the scaling properties of coarsening systems quenched to a critical phase and 
the results of numerical simulations of the clock model with p = 6 quenched into the KT phase. Sec. IVl contains the 
final observations and the conclusions. 



II. MODEL AND OBSERVABLES 



The p-state clock model is defined by the Hamiltonian 

H[<r] =~JJ2 3i ' B i = - J J2 cos & ~ i)> ( 3 ) 

<ij> <ij> 

where Si is a two-components unit vector spin pointing along one of the directions 9t = 2Trrii/p, with m S {1, 2, 

and < ij > denotes nearest neighbors sites i,j on a lattice. We will consider a square lattice in spatial dimension 

d = 2. This spin system is equivalent to the Ising model for p = 2 and to the XY model for p — ► oo. 

For p < 4 the clock model has a critical point separating a disordered from an ordered phase at T = Ti . For p > 5 
there exist two transition temperatures T\ and T^>T\ pj. For T <T\ the system is ferromagnetic, and for T > T2 
it is in a paramagnetic phase. Between these two temperatures, for T\ < T < T2, a KT phase Q exists where the 
correlation function behaves as G eq {r) ~| r | - ''( T ) with the anomalous dimension rj(T) continuously depending on 
the temperature. Both the transitions are of the KT type, namely the correlation length diverges exponentially as T\ 
or T2 are approached from the ferromagnetic or paramagnetic phase, respectively. Approximate analytic results [7j 
predict 

4tt 2 

T '/ J = W (4) 

and T2 to coincide with the KT temperature of the XY model T2/J = 0.95 (here and in the following we set the 
Boltzmann constant kg = 1). The exponent rj(T) is expected 7] to vary between ry(Ti) = 4/p 2 and ^(X^) = 1/4. 
Numerical simulations [9| are consistent with these predictions. 

A dynamics is introduced by randomly choosing a single spin and updating it with Metropolis transition rate 

w([a] -> [a'}) = min[l,cxp(-A£/T)] . (5) 

Here [a] and [a 1 ] are the spin configurations before and after the move, and AE = H[a'} — H[a\. 

We consider the protocol where the system is initially prepared in an high temperature uncorrelated state and then 
quenched, at time t = 0, to a final temperature Tf in the ferromagnetic phase or in the KT phase. The characteristic 
size L(t) grows until it becomes comparable with the system size and the new equilibrium state at Tf is globally 
attained. For an infinite system the final equilibrium state is never reached and L(t) keeps growing indefinitely. In 



3 



the late stage the power law {Q sets in. The characteristic length L(t) can be estimated from the knowledge of the 
two-points equal time correlation function 

G(r,t) = <<?,(*) ■ tj(t)) (6) 

where ai(t),aj(t) are spin variables at time t on two lattice sites whose distance is r, and (. . .) means an ensemble 
average, namely taken over different initial conditions and thermal histories. Due to homogeneity, G(r, t) does not 
depend separately on i an j but only on r. Enforcing this property we will numerically compute the correlation in 
the following as 

i je.J+ 

where i runs over all the N sites of the lattice and J+ is the set of four points reached moving a distance r from i 
along the horizontal or vertical directions. The methods to extract L(t) from G(r,t) depend on the scaling properties 
of G(r, t) and differ if quenches in the ferromagnetic or in the KT phase are considered, as discussed in Sees. IIIIIIVI 
The two time quantities that will be considered in this paper are the autocorrelation function 

C(t,s) = (a i (t)-a i (s)) (8) 

and the integrated (auto)response function, or zero field cooled susceptibility 



X(t,s)= I dt'R(t,t'). (9) 

The quantity 



(10) 

h=0 



a being a generic vector component, is the response function associated to the perturbation caused by an impulsive 
magnetic field hi switched on at time t' < t. Recently, new efficient methods for measuring the response function 
without applying the perturbation have been introduced |lOllllLll^ |. In the following we will use the one derived in |12| . 
For spin systems sub ject ed to a Markovian dynamics an out of equilibrium generalization of the fluctuation dissipation 
theorem was derived jl3j , relating the response functions to particular correlation functions of the unperturbed system. 
For the integrated response function © it reads 



T X (M) = 5 



C(t,t)-C(t,s)- / (ai(t) ■ Bi{t'))dt' 



(11) 



where 



B i [a] = -^ / (a i -a' i )w([a]^[a'}). (12) 



In this equation [a] and [a'] are two configurations differing only by the spin on site i, taking the values <7j and a\ 
respectively. The notation (<Xj(t) • Bi{t')) in Eq. (|TT]) means the average Ylw [a'] ^* ' ^i[ a ' / ]p([ a ']' *i [ <t ']>*')p([ <j/ L*')> 
where p([a'],t') is the probability to find the configuration [a 1 ] at time t' and p([a},t; [a'],t') is the joint probability 
between [a] at time t and [a 1 ] at time t' . 

Eq. I|ll|) allows to compute the integrated response function by measuring correlation functions on the unper- 
turbed system, avoiding the complications of the traditional methods where a perturbation is applied, and improving 
significantly the quality of the results . 

III. QUENCHES TO T f < Ti 
A. General scaling properties 



Let us start considering quenches from an high temperature disordered phase to the ferromagnetic region. In this 
case one observes the growth of compact domains separated by topological defects such as interfaces or vortices (see 
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FIG. 1: 

(Color online) C a g(t,s) obtained with the two methods described in the text is plotted against t — s for s = 
100, 200, 400, 800, 1600, 3200. 

Fig. [21 • As a consequence, a sharp distinction can be made in the late stage between spins belonging to the interior 
of domains from those pertaining to the defects. The interior of the domains very soon attains local equilibration 
in one of the broken symmetry equilibrium phases at Tf, whereas the degrees of freedom around defects are out of 
equilibrium and are responsible for the aging of the system. Observables such as C(t,s), x{^i s ) an d G{r,t) take an 
addictive structure 0] 

C(t,s)=C st (t-s)+C ag (t,s), (13) 

and similarly for \(ty s) and G(r,t). In the following, for clarity, we will focus on C(t,s), but similar considerations 
hold for the other quantities. The dynamics of the spins in the bulk of domains provide the equilibrium contribution 
C s t{t — s) while what is left over, C ag (t,s), accounts for the aging behavior. Since equilibrium dynamics is well 
understood, the behavior of C s t(t — s) is generally well known. In particular, at Tf = equilibrium dynamics is frozen 
and C s t(t — s) = 0. On the other hand, much interest is focused on the aging part of the aforementioned observables, 
which is less understood. This contribution can be isolated by subtracting C s t(t — s), computed in equilibrium, from 
C(t, s). However, for models with a discrete symmetry, it is computationally much more efficient to resort to a different 
method. This amounts to study a modified system where Tf in the transition rate JSJ is set equal to zero if the spin 
<7i to be updated belongs to the bulk, namely if it is aligned with all its neighbors. Since the bulk degrees of freedom, 
which alone contribute to C st (t — s), feel Tf — 0, and C st {t — s) = at Tf = 0, by computing observables with this 
modified dynamics one isolates the aging term leaving other properties of the dynamics unchanged |l5j | . In order to 
check this we have computed C(t,s) and C s t(t — s) with the Glauber dynamics in a system with p = 3 quenched 
to Tf = 1 < T\ and in an equilibrium system at the same temperature. In Fig. Q we compare C ag (t,s) obtained 
by subtraction of these quantities through Eq. (|13f) (symbols) and the same quantity obtained directly by means 
of a quench simulation with the modified dynamics (continuous lines). This figure shows an excellent agreement, 
confirming that the modified dynamics is efficient and accurate. In the remaining of this Section, therefore, we will 
always present results obtained with this modified dynamics. 

In the late stage of the evolution, after a characteristic time t sc when L(t) is much larger than all other microscopic 
lengths, dynamical scaling is obeyed Accordingly, for the correlation function one has 

G ag (r,t) = M 2 g(x), (14) 

where x — r/L(t) and M is the equilibrium magnetization at Tf. In systems with a discrete symmetry and sharp 
interfaces, a short distance behavior (x <C 1) of the type 1 — g{x) ~ x is found Q, ||J, namely a Porod's tail 
G(k,t) ~ fc~( d+1 ' in momentum space for large k (k ^> L{t)~ x ). This is known to be true, in particular, for the clock 
model, for all p < oo, although the whole form of g(x) depends on p In systems with a vector order parameter 
and an O(N) symmetry one has a generalization of the Porod's law O, G(k,t) ~ k-( d+N l 

From Eq. 114f) one can extract a quantity L\(t) proportional to the typical domain size L{t) from the condition 
g(x) = \, namely as the half-height width of G ag (r,t). Alternatively, a characteristic length L-^it) can be extracted 
as L 2 {t) = JdrG ag {r,t). Clearly, if Eq. (JUJ) holds, Li{t) oc L 2 (t). 

The size of domains can be related to the density of defects p(t) . For a coarsening system where topological defects 
are only interfaces, such as the Ising model, one has a power law behavior p(t) oc t~ s , with Jl| 

5 = l/z. (15) 

This result does not ap ply to systems with different defects, such as vortices or others. In the case of vector O(N) 
(with N > 2) model [l,ll7| one has 

S = 2/z, (16) 

and logarithmic corrections for N = 2. In these cases p(t) provides an indirect, alternative method for the deter- 
mination of L(t), and hence of z. For the clock model, where defects are interfaces and vortices, neither Eq. I|15|) 
or Eq. l|l(jD can be straightforwardly used. However, a simple inspection of the configurations (see Fig. |2Jl suggests 
that, since vortices are point like, their contribution to p(t) must be negligible in the late stage. Therefore we expect 
Eq. Q15[l to be obeyed asymptotically. 

The dynamical exponent z is believed to be universal for quenches to the ordered phase Tf < T\: the same value 
z = 2 as for the Ising model is expected for every value of p |5j,|6( in the clock model and for every N in O(N) models. 
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Coming to two-times quantities, the aging part of the autocorrelation function is expected Q, Q] to scale as 

C ag (t,s) = h(y), (17) 

with y — t/s and h(y) ~ y~ x / z for y ^> 1. The exponent A is believed to be the Fisher-Huse exponent which regulates 
the large t decay of the initial condition autocorrelation function C(t, 0) ~ t~ x / z . In the Ising model one has A = 5/4. 
We are not aware of a systematic study of this exponent in the XY model. In [H| it is argued that this exponent 
depends on Tf and, for the particular case Tf = 0.3 the value A = 0.54 is found. For the integrated response function 
scaling implies 



For p = 2 the scaling function behaves as 



Xa g {t,s) = s- a *f{y). 



f(y) ~ y~ 



(18) 



(19) 



for y ^> 1. Regarding the exponent a x , analytical calculations in solvable scalar models or in the large- N model 19, 2 
find the following dependence on dimensionality 



for 



d < d, 



u 



dy = s 5 with log corrections for 
for d > du, 



(20) 



where du is the lower critical dimension of static critical phenomena and du is an upper dimension that turns out to 
be du — 3 or du — 4 for systems with a discrete or continuous symmetry. This expression shows that the response of 
coarsening systems depends on dimensionality in a non trivial way. Numerical simulations 

El El E2 of scalar and 

O(N) vectorial systems, with conserved and non conserved order parameter, are consistent with Eq. (|20|l . The value 
of a x has never been investigated for systems with a discrete symmetry and a degeneracy of the ground state larger 
than p — 2, as in the case considered in this paper. 

Notice that, for a given dimensionality d < du, Eq. I|20|) predicts a different exponent for systems with a continuous 
or a discrete symmetry. In the case d = 2, for instance, for the Ising model Eq. (1201) gives a x = 1/4 while for the XY 
model a x = 0. Therefore, for the model under investigation, the value of a x may be used to discriminate between the 
Ising and XY non-equilibrium universality classes. 

Finally, let us recall that the scaling behaviors (|14I17I18|) are only expected asymptotically. Since numerical simula- 
tions can only access a finite time region, preasymptotic effects may be present. In particular in numerical simulations 
of the Ising model with a non-conserved order parameter, one usually observes an effective exponent l/z e ft ~ 0.48 
in place of 1/z = 0.5. The integrated response function has been also shown to be affected by corrections to 

scaling. These can be conveniently discussed in terms of the effective exponent, defined as 



a eff(V, 



dlllXa,g(t, s) 



d In , 



(21) 



With a scaling form such as l|18|l . one would have a e ff(y, s) = a x , independent of y and s. However, if preasymptotic 
effects are present, the effective exponent takes a value which depends both on y and s. For p = 2 it was shown [l5l | 
that, because of this, a e ff(y, s) is found in numerical simulations in the range 0.25 < a e ff(y, s) < 0.28. 



B. Numerical results 



In the following we will present the numerical results. Setting J = 1, for each case considered we simulated a 
square lattice of size 1000 2 with periodic boundary conditions and an average over 100 realizations was performed. 
Statistical errors, when not explicitly plotted in the figures, are comparable to the thickness of the lines. 

For p = 3 there is a ferro-paramagnetic transition at T\ ~ 1.326 while for p = 6, according to Eq. |0J, one has 
Ti ~ 0.645. We performed a series of simulations of quenches to Tf < T\, with Tf = 1/2 or Tf — 1 for p = 3 and 
p = 6, respectively. Typical configurations of domains in the late stage are shown in Fig. El Notice the simultaneous 
presence of interfaces and vortices. These are defined analogously to those of O(N) models: on encircling a vortex 
the order parameter rotates by ±2ir (although in the clock model rotations are obtained by discrete steps). While 
for p — 3 vortices and interfaces between different phases are all energetically equivalent, in the case p = 6 one has 
the additional feature of different kinds of vortices and interfaces. Let us consider a domain characterized by having 
all the spins pointing along the direction 8 = 2itn/p. Spins belonging to the domain are characterized by having the 
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FIG. 2: (Color online) Configuration of the system with p = 3 (left) and p — 6 (right), at t = 3200. 



same value of rii, rij = n. This domain can be separated by an interface from another domain where spins point along 
a different direction 9 = 2irm/p. Clearly, interfaces between domains of contiguous phases, namely with n = m ± 1 
are energetically less expensive then the others. The more energetically expensive interfaces are eliminated faster 
from the system (they are already practically absent in Fig. [2}. This fact influences considerably the topology of the 
growing structure. For instance, one clearly observes that between two domains of non contiguous phases, say with 
n = m and n = m + 2, a thin slab of phase n = m + 1 is interposed in order to minimize energy. 

Analogously one observes also different kinds of vortices. Points where six phases meet are energetically favored, 
but vortices where a lower number of phases meet can also be present. Moreover there are also points where four 
(or more) domains meet, but two of them belong to the same phase: Encircling such points one may enter domains 
characterized respectively by, say, the sequence n, n + 1, n + 2, n + 1 again. Clearly, encircling the most energetically 
favored vortices one finds all the phases according to the sequence n, n+ 1, n + 2, n + 3, n + 4, n + 5, n + 6 (or in reverse 
order). As for the interfaces, the high energy vortices are quickly removed and a typical late stage configuration, as 
that of Fig. |2 contains practically only the lowest energy vortices. The presence of all these different kind of defects 
in the system is possibly the origin of the long lasting preasymptotic effects discussed below. 
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FIG. 3: (Color online) Comparison between Li(t),L<2(t) and p(t) 1 for p = 3 (left) and p = 6 (right). 



FIG. 4: (Color online) Data collapse of G ag (r,t) against x = r/Li(t) for several times t„ generated from t„ =Int[exp(n/2) + 1] 
with n ranging from 13 to 22 and p — 3 (left) or p = 6 (right). 

A comparison between L\{t) , L,2(t) and is shown in Fig. EI After an early stage when domains are formed 

and scaling does not hold, Li(t),L2(t) and p(t) _1 start growing with an approximate power law behavior and for 
long times one has L\(t) oc L 2 {t) oc /5 -1 (t) (for p — 6, /o(t) _1 does not obey a power law behavior in the range of 
simulated times. However, for the longest simulated times the effective exponent seems to approach a value roughly 
comparable with that of L\[t) and L 2 (t)). This implies that Eq. I|15fl is obeyed asymptotically. We recall that for a 
system containing only vortices, such as the XY model, one would instead expect the relation Regarding the 

coarsening exponent, in the decade 10 4 — 10 5 for p — 3 we measure l/z e ff — 0.486 ± 0.002, l/z e // = 0.484 ± 0.002 
and l/z e ff = 0.478 ± 0.002 from L x (t), L 2 {t) and /o(t) -1 , while for p = 6 we get l/z eff = 0.467 ± 0.003, l/z eff = 
0.474 ± 0.003, 1/Zeff = 0.450 ± 0.008. These values (apart from the last one which is evidently a still preasymptotic 
effective exponent) are compatible with the value 1/z — 1/2 of the Ising model. The different initial behaviors of 
L\{t), L 2 (t) and p(i) _1 , signal that preasymptotic effects are present up to very long times. This is probably related 
to the presence of different types of defects. Note also that, at the longest time considered, the density of defects with 
p = 6 is more than three times larger than with p = 3. 

In FigQJwe test the scaling form (|14|) of the equal time correlation function. We plot G ag (r,t)/M 2 against x — 
r/Li(t) for several values of t in the two decades range [6.4 • 10 2 — 6.4 • 10 4 ]. The data show a good collapse on a single 
master-curve g(x). For small x the Porod's behavior 1 — g(x) ~ x is very neatly observed. 

We turn now to consider two time quantities. In Fig[S]the autocorrelation function is plotted against y for different 
values of s in the range [100 — 3200]. As already observed when discussing Fig. |3 scaling is only approximatively 
obeyed in this regime. This is mirrored by C ag (t, s): Full data collapse is not found. Regarding Fig. El we also noticed 
that scaling improves as time gets larger and it is reasonably well obeyed for the longest simulated times. The same 
conclusion can be drawn, for p = 3, from C ag (t, s). Indeed, in Fig[5]one can observe that the data collapse improves 
pushing s and y to larger values. In fact, although the data collapse of the curves is poor for small y it gets better 
increasing y and, for y > 10, all the curves collapse. Moreover, the quality of the collapse improves also increasing 
s. Indeed, while the curves with small s do not collapse (except, as anticipated, for y as large as y > 10) the two 
curves with the largest values of s (s = 1600 and 3200) practically coincide for all y > 2. For p — 6 the collapse 
is worse. Coming to the asymptotic behavior of the scaling function h(y) ~ y~ x l z , it is numerically too demanding 
to reach the asymptotic large-y region with the values of s considered in Fig. [S] Then, we have evaluated A from 
the large t behavior of C(t,0). This quantity is shown in the insets of Fig. [5] In the range t € [4 • 10 4 , 10 5 ] we find 
X/z = 0.61 ± 0.01 and X/z — 0.57 ± 0.01 for p — 3 and p = 6. For p = 3 this value is consistent with the value 
X/z = 5/8 = 0.625 of the Ising model, keeping also into account that the effective exponent we measure is still slightly 
increasing at the longest simulated times. For p = 6 the measured exponent is somewhat smaller than for p — 3 and 
for the Ising model, but the fact that it keeps still growing at the longest simulated times suggests that asymptotically 
the same value A = 5/8 = 0.625 could be obtained. This results suggest that there could be a unique non-equilibrium 
universality class for every 2 < p < oo. This hypothesis can be tested by considering the integrated response function. 
This quantity, for p = 3, is plotted against t in Fig. [21 showing a marked dependence on s. Starting from zero at t = s, 
Xag(t, s) reaches a maximum at t ~ 2s and then decreases to zero with a power law behavior. A similar behavior is 
observed for p — 6. According to Eq. (|18f) . by fixing y to a certain value and varying s or equivalently t, the data 
should follow a power law with exponent — a x . As an example, the points corresponding to y = 4, which have been 
marked with stars in the log-log plot of Fig. ©, are approximatively aligned on a straight line of slope 0.26 ± 0.01. 
A similar analysis can be performed for every value of y. 



FIG. 5: (Color online) C(t, s) is plotted against y for p = 3 (left) with s = 100, 200, 400, 800, 1600, 3200 and for p = 6 (right) 
with s = 400, 800, 1600, 3200. In the insets C(t, 0) is plotted against t. 
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FIG. 6: (Color online) TfXag{t, s) is plotted against t for fixed values of s. The dashed line is the power law t 



FIG. 7: (Color online) Tf\ a g{t, s) is plotted against s for p = 3 (left) and p = 6 (right), with fixed values of y (y = 
2, 4, 6, 8, 10, 12, 14, 16, 18 from top to bottom). Numerical values are marked with error bars, continuous lines are guides for the 
eye. 



In order to make a quantitative analysis of this exponent, and to detect preasymptotic effects, in Fig. (JZJ we plot 
Xag(t, s) for fixed values of y against s. According to Eq. j^Jl the slope of these lines is a e //(y, s) and, if scaling ljT%|) 
holds, one should find a e ff(y, s) = a x . Preasymptotic effects, instead, introduce a weak dependence of this exponent on 
s and y. Apart from the curve y = 2, which corresponds to very early times, the slopes of all the curves are compatible 
with Eq. JSUJl, namely with a x = 1/4 (for p = 3 in the range s G [800 — 3200] we find a e ff(y, s) — 0.25 — 0.30, depending 
on y. For p — 6 the effect of preasymptoticity is smaller and one finds a e //(y, s) < 0.27 for every value of y and s). 
This pattern of behavior of a e //(y, s) is similar to what is observed in the Ising model where a e //(y, s) ranges in the 
interval, a eff (y, s) = [0.25 - 0.28]. 

The data collapse of s^Xagl^, s) vs y expected from Eq. i|18[l is shown in Fig. EI F° r P — 6 the collapse of the two 
curves with the largest s (s = 1600,3200), is good at sufficiently large y (y > 4). It is poorer for p = 3. Note also 
that the asymptotic behavior f(y) ~ y~ 1//4 for y — > oo is well obeyed, consistently with Eq. H19fl . again confirming 
a x — 1/4 and ruling out the value a x = appropriate to the XY model. These results strengthen the conclusion that 
the clock model below T\ does not belong to the XY universality class and that the non equilibrium universality class 
is the same for all 2 < p < oo. 



Let us consider the model with p = 6 quenched to the critical region T\ < Tf < Ti. In the late stage the correlation 
function obeys Eq. J5J. From this equation one has 



IV. p = 6: QUENCHES TO Ti < Tf < T 2 




(22) 



L(t) can then be extracted as 



L{t) oc I{t) 



(23) 



The autocorrelation function obeys 0, |2^| 



C{t,s) = (t-s + to) 




(24) 



FIG. 8: (Color online) TfS 1 ^ 4 Xag(t, s) is plotted against y for fixed values of s for p = 3 (left) and p = 6 (right). 
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where to is a microscopic time. Neglecting to for t — s ^> to, C(t, s) can be rewritten in scaling form 

{d-2 + r,) 



C(t,s) 



-h(y) 



(25) 



where h(y) = (y — l)~( d ~ 2+,1 ^ z h(y), with the property h(y) ~ y~ x / z for y ^> 1. Notice that, when using the y variable 
in the large-s limit, the condition t — s 3> to for the validity of Eq. i|25[l becomes y > 1. For all y > 1, then, one should 
expect data collapse of the curves s^ d ~ 2+rt ^ z C{t, s) against y for different choices of s. 
The response function is given by 0, |23| 



(26) 



(d-3 + 7}) - ft 1 

R(t,t , ) = {t-t , +to)--^ M fr 1 



Splitting the integral of Eq. Q into two integration domains, introducing the arbitrary number e, the integrated 
response function can be written as 



X(t,*)~i- 



du (1 — it H ) » J\ u )+ I du(l — u-\ ) » /(u) 

5 E 



(27) 



(28) 



where u = t! ft. For large times t one can choose t /t «C e <C 1/2. The condition £ /£ <C e allows one to neglect to/t 
with respect to 1 — u in the first integral, whose value can then be evaluated as Ii(y, to/t, e) ~ F(l — e, to/t) — F(y, 0), 



where dF(u,t /t)/du = (l-u + 



f{u). Let us consider now the second integral l2(to/t,e) — F(l,to/t) — 



F(l-e,t /t). Here, since e < 1/2, one can set /(«) ~ /(l), so that F{l,t /t) ~ (t /t)-( d - 2+, ')/ 2: /(l)z/(d - 2 + r?). 
One then arrives at 



(d-2+77) 
X(i,s)=t — /(y) 



/(I) .- 



(29) 



where /(y) = —F(y,0). Letting t — > 00, x(t, s) must converge to the equilibrium susceptibility whose value is given 
by the fluctuation dissipation theorem, \ eq = T^T . This leads to the identification of the last term on the r.h.s. of 
Eq. with Xeq , and Eq. can be cast as 

X(t,s)- X e g ~s- a *f(y), (30) 

where a x = (d — 2 + r\)/ z. Notice that, differently from the case of quenches in the ordered phase, here the exponent 
a x is directly related to the equilibrium critical exponents rj and z. We recall that, in the KT phase, the exponents 
77, z, A depend on temperature. 

It is interesting to discuss the parametric plot of \(t, s) against C(t, s). Since C(t, s) is a monotonically decreasing 
function of t, this time can be re-parametrized in terms of C, obtaining %(i, s) = x(C, s). This quantity is important 
because, if appropriate conditions are satisfied, its large-s limit 



X (C) = lim X (C, s) 

S — M30 

provides a connection between static and dynamic properties through the relation 

d 2 x(C) 



f 



dC 2 



(31) 



(32) 



C=q 



where P{q) is the overlap probability function of the equilibrium state at T — Tf. As discussed in |14| . a universal 
linear relation 



TfX(C) = T fX 



eg 



c 



(33) 



as for equilibrated systems, is expected for quenches to a critical point or into a critical phase, although the system 
is aging for any finite time. 
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FIG. 9: (Color online) Configuration of the system at t = 3200 after a quench to Tf = 0.76. 



FIG. 10: (Color online) The behavior of I(t). The dashed line is the power law t 1 



A. Numerical results 

In this section we present results of simulations of the model with p = 6 quenched to T\ < Tf = 0.76 < T^. A 
typical configuration of domains in the late stage is shown in Fig. [!|] In this case, there are no compact domains. 

The quantity lit) is shown in Fig. ^| Here one observes that the power law behavior sets in very early. This 
implies, through Eq. I)22[l. that also Lit) has a power law growth. The effective exponent has a small tendency to 
increase as t gets larger: We measure an exponent 0.34 in the decade 10 2 — 10 3 and 0.35 for t > 10 4 . The exponents r\ 
and z are known numerically 9]. At Tf — 0.76 their measure gives rj = 0.17 and z = 2.18, yielding (3—d—r))/z = 0.38. 
This number is consistent with the value 0.35 obtained from I(t) by means of Eqs. (|23|l . taking also into account that 
the effective exponent we measure is still increasing at the longest simulated times. 

In Fig. ^2 we test the scaling form of the equal time correlation function. We plot r v G a g(r, t) against x for 
several values of t, where L(t) is computed through Eq. i|23|l . The data show a very good collapse on a single master- 
curve g(x). Notice that, as expected, Porod's behavior at small x is not observed, due to the non compact nature of 
the domains. 

We turn now to consider two time quantities. In Fig. 1121 the the autocorrelation function is shown. There is a 
tendency to a better data collapse for larger times, implying that the scaling symmetry is still not exactly obeyed. 
For the two largest values s = 50 and s = 100, however, the collapse is rather good. 

Let us consider the integrated response function, that is shown in Fig. 1131 Here one observes an analogous situation: 
the collapse expected on the basis of Eq. (|3"U|l is rather good for the two largest values of s. 

In Fig. [21 the parametric plot of x(C, s) is shown. For the largest values of C, x(C, s) obeys Eq. l|3"3"|) . As C is 
decreased the curves flatten and Tfx(C,s) lies below the asymptotic curve l|33[) . However, in the limit of infinite 
times t — > oo , which corresponds to C —> 0, each curve must necessarily obey Eq. I|33l) . since x(t, s) approaches the 
equilibrium value Xeq — Then, moving toward C = 0, at some point the curves become steeper in order to 

meet the value Xeq at C = 0. Changing s, the same qualitative behavior is observed, but the curve gets higher, slowly 
approaching the asymptotic form (|33|) for all values of C in the large s limit. This pattern of %(C, s) is analogous 
to what observed in the spherical model quenched at the critical point and is expected in full generality whenever a 
system is quenched to a critical point or to a critical phase . It must be noticed that the convergence to the trivial 
form (|33|l is very slow because it is regulated by the rather small exponent r)(T) 0]. Since the exponent r?(T) at 
the lower transition temperature is expected to behave as r](Ti) ~ 4/p 2 , the asymptotic behavior can be arbitrarily 
delayed increasing p. The simulations presented in this section have been planned out in order to show at least a 
glimpse of this convergence. As discussed in Ref. [yji previous studies of the KT phase of the XY model, for which 
the same asymptotic form (|33|l is expected, interpreted a preasymptotic non trivial form of y[C, s) analogous to the 
one of Fig. E| as a reminiscence of the parametric plot of the d — 3 Edwards- Anderson model jl8| or used it to infer 



FIG. 11: (Color online) Data collapse of r v G(r, t) against x for several times (t = 666, 1097, 1809, 2981, 4915, 8104, 13360, 17155) 
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FIG. 12: (Color online) Data collapse of C(t, s), for s = 10, 20, 50, 100 (from bottom to top). 



FIG. 13: (Color online) Data collapse of x(t> s ) m the quench at Tf = 0.76, for s = 10, 20, 50, 100 (bottom up). 



the asymptotic value of the fluctuation-dissipation ratio |27j . The simulations presented here clearly show, instead, 
that this pattern is preasymptotic and the data are consistent with a convergence to the expected trivial limiting 
form 

V. CONCLUSIONS 

In this paper a rather general numerical investigation of the off-equilibrium dynamics of the clock model after a 
temperature quench has been carried out. We have considered both quenches into the ordered phase Tf < T\ for 
systems with p = 3 and p = 6 and a quench to the critical, Kosterlitz-Thouless phase T\ < Tf < T%, for p — 6. In 
all these cases we analyzed the behavior of one-time quantities, such as the equal time correlation function or the 
characteristic length, and two-times quantities, such as the integrated response and the autocorrelation function. This 
study provides a quite general scenario of the scaling properties of the dynamics and allows the comparison with the 
behavior of other well studied coarsening systems such as the Ising model or the XY model, which correspond to the 
case p = 2 and p = oo, respectively. We find that dynamical scaling is obeyed in all the cases considered. In the 
ordered region, the dynamical exponents are the same of those of the Ising model. While the result z = 2 for the clock 
model was already well known |5J, 6j the values of the exponents A and a x have been measured for the first time and 
deserve some considerations. These values suggest that the clock model belongs to the non-equilibrium universality 
class of the Ising model. Moreover, finding the same values of the exponents for both p — 3 and p = 6 implies that the 
system is in the same equilibrium universality class for all p < oo. Finally, the value a x = 1/4 fits with the general 
phcnomenological formula (|20|l for coarsening systems. This strengthens the idea that the non trivial dimensionality 
dependence of a x predicted by Eq. (|20H . may have a general validity for coarsening dynamics. 

It has been recently proposed |22j,|28| that, limited to the case of systems where topological defects are exclusively 
domain walls, Eq. (|20|l may be related to the dynamical roughening of the interfaces. The fact that the value of a x is 
the same also in the clock model implies that the asymptotic contribution of vortices to the response function is not 
important, at least at the level of the exponent a x , and suggests that the behavior of the interfaces is the unifying 
feature that makes systems with different degeneracy p fall into the same universality class. 



FIG. 14: (Color online) Parametric plot of x(*. s ) vs C{t,s) in the quench at T f = 0.76, for s = 10,20,50,100 (bottom up). 
The dashed line is the asymptotic curve Tf\(C) = TfXeq — C, Eq. 1331 . 
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